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I. INTRODUCTION 

In this paper we demonstrate tight binding (TB) mod- 
els for iron with and without interstitial hydrogen im- 
purites at the concentrated and dilute limits. Although 
there is a large number of existing classical potentials 
which are certainly of great importance and usefulness, 
they all suffer from a particular drawback in that the un- 
derlying classical EAM-type models for pure Fe, with one 
apparent exception^ fail to predict the known core struc- 
tures of screw dislocations*^ On the other hand tight 
binding models abstracted into bond order potentials cor- 
rectly predict core structures in agreement with first prin- 
ciples calculations. 2 Ultimately one of the many goals is 
to study how interstitials form atmospheres around dislo- 
cations and impede or enhance flow through mechanisms 
such as hydrogen enhanced local plasticity 3 (HELP) and 
so it is essential that dislocation core structures are cor- 
rectly predicted. A further slightly disturbing feature of 
the classical models is the truly vast number of parame- 
ters involved which have to be fitted to a very large train- 
ing set of data. Here in common with the approach to 
classical model fitting, 1 we first construct models for pure 
iron and then go on to make models for hydrides with- 
out further adjusting the Fe-Fe interaction parameters. 
But in contrast we try to keep the number of parame- 
ters and fitting targets to a minimum and focus on the 
ability of the models to predict those properties that are 
normally included in the training sets in the construction 
of classical potentials. We would argue that this is pos- 
sible because the TB approximation comprises a correct 
quantum mechanical description of both magnetism and 
the metallic and covalent bond and so the correct physics 
is built in from the start. That being the case, we do not 
expect the theory to be over sensitive to the choice of pa- 
rameters and indeed in the procedure we describe below 
a large number of equally useful models is thrown up. 

The structure of this paper is as follows. In section [TT1 



we revisit the tight binding approximation and discuss its 
parameters and their environment dependence, or screen- 
ing. We describe two models for pure Fe in section IIIII 
which are fitted to properties of bulk bec a-Fe and used 
to predict properties of fee 7-Fe and hep e-Fe, as well 
as surface and vacancy formation energies in a-Fe. In 
section IIVI we augment one of these models with Fe-H 
interactions which we fit to the properties of four mono- 
hydride FeH phases, and test against adiabatic potential 
surfaces. We then proceed to the dilute limit of H in Fe in 
section fVl without further adjustment of parameters and 
use our model to predict segregation energies of H to in- 
terstitial sites, vacancies and surfaces of a-Fe. By and 
large, we find remarkable agreement with published ex- 
perimental results and ab initio calculations. We discuss 
our models and conclude in section IVll 

II. THE TIGHT BINDING APPROXIMATION 
AND TRANSFERABILITY 

A. Distance scaling and range of the hopping 
integrals 

There is no need to rehearse the tight binding ap- 
proximation in any detail here. Recently Paxton and 
Finnis 4 constructed tight binding models for magnetic 
Fe and Fe-Cr alloys and details can be found there as 
well as in many other publications. 5-9 However we do 
wish to make some preliminary remarks. The scheme 
that we use is the self consistent Stoner model for itiner- 
ant ferromagnetism 8 and goes beyond the fixed moment 
and rigid band approximations. The connection between 
tight binding theory and the first principles local spin 
density approximation (LSDA) to density functional the- 
ory (DFT) is now well established] 7 ' 10 ' 11 TB is compu- 
tationally several orders of magnitude faster than LSDA 
because the hamiltonian is constructed from a look-up 
table of parameterized hopping integrals, h, and possibly 
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FIG. 1. (color online) Energy bands for bee Fe at its experimental lattice constant, 2.87 A. The coloring is such that s character 
is green and d character is blue. The Fermi energy is indicated by a horizontal line. The upper panels are majority and the 
lower panels minority spin bands. Far left is the tight binding d-band model and in the center our non orthogonal sd model. 
To the right are bands calculated in the LSDA-GGA. 
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overlap integrals, s. These are conventionally written 
in Slater and Koster's notation^ as ssa, sda, dda, ddn, 
ddS. Central to a tight binding model is the way in which 
these integrals scale with bond length. In this work we 
will use^i2ri& 

h(r) = h e-« r (1) 

and similarly for overlap integrals, when used, 

s(r) = s e- qr . (2) 

The alternative is to use the power law scaling, h ~ r~ n , 
demanded by canonical band theory^ - — There is no 
strong argument to prefer one over the other; in fact 
by equating^ logarithmic derivates of h{r) at, say first 
neighbors at a distance ro , we have n = qrg and in the bee 
structure of Fe q « 1 a.u. corresponds to the canonical 
n — 5 (see table HI below). 

This brings us to a well known paradox of tight binding 
modeling namely that the decay of the hopping integrals 
is known a priori from band theory, which may render 
them longer ranged than is desirable. A well known ex- 
ample is the group IV semiconductors where by analogy 



with the free electron bands, to reproduce the volume 
dependence of the bandwidth the hopping integrals must 
scale with n = 2i2i This scaling is bound to lead to very 
long ranged hopping integrals; on the other hand it is 
known that the first neighbor approximation is the right 
one, and attempts to include further neighbors fail^ For 
many purposes it is adequate simply to cut off the in- 
teractions between first and second neighbors, but this 
can lead to difficulties in work on complex defects or in 
molecular dynamics. An elegant solution was provided 
by Goodwin, Skinner, and Pettifor— (GSP) which cuts 
off a power law exponentially beyond some chosen cut-off 
distance, r c . There are two drawbacks to this, (i) An ex- 
ponential decay can still lead to discontinuities in molecu- 
lar dynamics (as one still needs to impose a cut-off in the 
neighbor lists), (ii) The GSP form maintains the value 
but not the slope of the underlying power law at first 
neighbors. Therefore our preference is to retain the power 
or exponential scaling given by the canonical band theory 
and to choose two distances, r% and r c , between which 
to smoothly augment the interaction to zero. This can 
be achieved by matching value, slope and curvature at r*i 
and at r c with a fifth degree polynomial which replaces 
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the hopping integral in that ranged We show our hop- 
ping integrals thus augmented at figure [5] in section [TV] 
below, where we discuss this matter further. 

B. The pair potential, transferability, and non 
orthogonality 

The hopping integrals provide an attractive force, 
which in the conventional tight binding models is bal- 
anced by a repulsive pair potential, which here may take 
the form 

0(r) = B ie - pir - B 2 e~ P2r (3) 

in which, as suggested by Liu et al&, both Bi and B 2 
are positive. This potential is expected to be repulsive 
at short range but is not positive for all r (see Fig. [4j 
below). 

An additional non pairwise repulsion is provided if it 
is chosen to make the model basis non orthogonal. This 
may give a number of advantages^ One is, that it is 
widely believed that non orthogonality confers a greater 
transferability to the model^ By this is meant that a 
model constructed for a particular crystal structure is 
less likely to fail when transferred into a situation of dif- 
ferent crystal structure or increased or reduced coordina- 
tion. We will wish to focus critically on this aspect of our 
models below. It is instructive at this stage to recall that 
by its very construction the tight binding approxima- 
tion discards all three center terms in the hamiltonian. 9 
On the one hand the canonical band theory shows that 
these, like non orthogonality, are of second order in the 
band widthi 19 i 26 ' 27 On the other hand Tang et al~ and 
Haas et al~& took the important step of proposing en- 
vironment dependent hopping integrals. In this empir- 
ical scheme the hopping integral between two atoms is 
modified in the close proximity of a third atom — in the 
extreme limit this third atom may approach the two cen- 
ter bond, generally speaking weakening or "screening" it, 
and eventually come in between the two atoms. Whereas 
the screening was first described by an empirical for- 
mula, Pettifor succeeded in deriving the Tang et al^ 
expression from the Lowdin transformation of the non 
orthogonal hamiltonian.— In particular he showed that 
sd overlap matrix elements in pure transition metals pro- 
vide this "screening" of the two center bond. Therefore 
rather than adopting explicit environment dependence 
as is done in recent bond order potentials; 31 ^ 32 we retain 
the two center approximation and employ non orthogonal 
models to account for the screening. 

C. The choice of parameters 

A related and highly significant finding 30 is that hop- 
ping integrals extracted from an LSDA hamiltonian cal- 
culated using the tight binding LMTO-ASA method 27 
are discontinuous between first and second neighbors in 



bec transition metals. These discontinuities are described 
consequently by the screening — a feature of the geometry 
of the bec lattice — leading to the analytic form of Tang et 
al££- and Haas et aZ,™ The point we wish to raise here 
is that it became clear— that transferable hopping inte- 
grals may be extracted from an LSDA hamiltonian thus 
avoiding the usual need for fittingi^i^ . Of course there is 
no unique tight binding model for a given element since 
the LSDA hamiltonian is basis-set dependent. We do not 
adopt this approach here for two reasons. First, the hop- 
ping integrals deduced from the LMTO-ASA 9 -^ derive 
from a hamiltonian whose on-site matrix elements are 
strongly volume dependent whereas in the tight bind- 
ing approximation these terms are volume independent 
and hence any volume dependence of the electronic struc- 
ture must be taken up by the scaling law (|TJ) . Second, if 
the hopping integrals and their scaling are taken from ab 
initio bandstructures without permitting further adjust- 
ment, then essential properties such as elastic constants, 
lattice constants and structural energy differences may 
have to rely on the choice of pair potential placing a large 
burden on that part of the model which is the most ad 
hoc. 



III. MODELS OF PURE IRON 



A. Orthogonal d, and non orthogonal sd models 

Construction of a tight binding model for transition 
metals is quite straight forward if it is not required to 
take the parameters from first principles bandstructure 
calculations A Given that the scaling should be close to 
canonical, as should the ratios^ 

dda : ddir : ddS « —6 : 4 : —1 

it is simple enough to guess a set of hamiltonian and 
overlap matrix elements and adjust these until the re- 
sulting energy bands match reasonably closely those from 
the LSDA. In fact, in all that follows we have used the 
LSDA with a generalized gradient correction (GGA) of 
Perdew et alM- With the exception of data taken from 
the literature all our LSDA-GGA results are calculated 
using the full potential LMTO method^ Energy bands 
calculated in this way are shown on the right in Fig.[TJ A 
simple canonical d-band model produces the bands shown 
to the left of Fig. [T] In addition to the integrals already 
discussed, we require a Stoner parameter, /, which rep- 
resents an on-site Coulomb integral^ to achieve a split- 
ting of the up and down spins Furthermore since the 
canonical model omits the s-band which is occupied by 
roughly one electron^ it is necessary to fix a number 
of d electrons, Nd^^- This is both the most simple and 
most reliable model for transition metals*^ 3 ^— Never- 
theless for the present purposes we wish to extend this 
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TABLE I. Parameters of our tight binding models for pure 
Fe. The {h} and {s} are the ho and so of equations ([TJ 
and ([2]). All quantities are given in atomic Rydberg units 
(1 bohr = 0.529 A, 1 Ry = 13.61 eV). Note that in the min- 
imal basis, orthogonal d model, the number of d-electrons, 
Nd is a parameter.^ Both hopping integrals and pair poten- 
tials are smoothly cut off between distances n and r c . These 
are shown in units of the bcc lattice constant, a. Both pair 
potentials are cut off with n = 1.1 and r c = 1.4, that is, 
between second and third neighbors of the bcc lattice (see 
Fig. [4] below) . By expressing n and r c in units of a we im- 
ply that these scale with the lattice constant, for example in 
the calculation of the bulk modulus, so that in a perfect lat- 
tice first and second neighbors always see the "proper" pair 
potential (|3]) and hopping integrals. This also applies below 
(Fig. [|)J to energy-volume curves in FeH, both for first and 
second bcc neighbors and first fee neighbours. 
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model. In the interests of transferability and to account 
for the bond screening without explicit environment de- 
pendent bond integrals, we explore here the addition of 
an s-orbital to the basis, including sdo~ and ssa non or- 
thogonality. We will also give arguments for this ne- 
cessity when we come to the iron hydrides below. The 
resulting bands, again obtained by simple comparison by 
eye with the LSDA-GGA bands are shown in Fig.Q] Den- 
sities of states associated with the LSDA-GGA and tight 
binding models are shown in Fig. [5J 



Having obtained two sets of bond integrals, we pro- 
ceed to find parameters of the pair potential and we do 
this by adjusting the four parameters in ([3]) to the lattice 
constant and the three elastic constants of bcc Fe. This 
cannot be done exactly because of the restricted form 
of the pair potential. The parameter sets are given in 
table U and resulting properties are in good agreement 
with experiment or LSDA-GGA calculations as can be 
seen in table [TTJ Calculations, written in italics in ta- 
ble [TT1 have been done using the full potential LMTO 



FIG. 2. Densities of states of pure bcc Fe in the orthogonal d 
(a) and non orthogonal sd (b) tight binding models compared 
to the LSDA-GGA (c). Majority and minority spin states are 
shown in the upper and lower panels of each. The zero of 
energy is shifted to the Fermi energy. 
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method^ and elastic constants are all obtained at the 
theoretical lattice constants. Our calculated lattice and 
elastic constants are in general agreement with previous 
work 42 



B. Predictions of the models 

1. Magnetic moment and structural magnetic energy 
differences 

It should be noted that the two models we have de- 
scribed are rather intuitively obtained and so, apart from 
the pair potential it cannot be said that these are "fit- 
ted" in the sense of a classical potential. Hence the 
properties shown in table |H] are in essence predictions of 
the model, validating the underlying correctness of the 
tight binding theory. These predictions can be discussed 
in more detail by reference to Fig. [3] which shows the 
structural energy-volume relation in bcc and hep Fe bro- 
ken down into bandstructure energy and magnetic energy 
contributions^ Both models reproduce the essential fea- 
tures which are, (i) the rapid collapse of the hep magnetic 
moment under pressure; (ii) the slow decline of the bcc 
moment and (m) the stabilization of bcc over hep being 
a result of the magnetism. The role of the pair poten- 
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FIG. 3. (color online) Contributions to the energy-volume relation in the orthogonal d (left) and non orthogonal sd (right) 
tight binding models. The lower panel shows the volume dependence of the magnetic moment. The dotted lines refer to the 
hep crystal structure and show the rapid collapse of the moment under pressure. The solid line is the bec moment and may 
be compared with the circles which are LSDA-GGA calculations. The upper panel shows the bandstructure energy (blue) and 
the magnetic energy (green) and their sum in red. Solid lines refer to bec and dotted lines to hep. The pair potential energy 
favors bec in both models. A vertical line indicates the equilibrium volume in bec Fe. 
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tial warrants explanation here. Fig. 0] shows equation (|3]) 
plotted using the parameters of our two models from ta- 
ble HI It might well be supposed that the stability of 
the bec phase compared to the hep is an artefact of the 
negative region of </>(r) falling at the second neighbors 
of the bec structure, while the 12 hep nearest neighbor 
distances fall in a positive region. This would be a valid 
criticism of our and Liu et al.'s~ models but is mislead- 
ing. In fact we find that we can easily make models 
that stabilize bec employing a pair potential that is pos- 
itive everywhere. In addition the stabilisation of the bec 
structure can be amplified by choosing larger Stoner I pa- 
rameter. We allow a larger moment in our orthogonal d 
model since it is known that the magnetic moment in 
bec Fe would be closer to 2.6/xb in the absence of sd and 
pd hybridization^^ Therefore the LSDA-GGA bec-hep 
energy difference is better rendered in that model (ta- 
ble |lll whereas we have chosen a value of / in our non 
orthogonal sd model that strikes a compromise between 
a smaller bec-hep energy difference having the benefit 
of a magnetic moment closer to the observed value. The 
real benefit of the form ^ is that it enables a sufficiently 
large value of the elastic constant C which otherwise ap- 
pears too small. It is well known that C can become very 
soft in bec metals and the values we obtain in table HJ are 
the best we can achieve after many trials with the other 
parameters and scalings in the models. Indeed in the 
model of Liu et al&, C is significantly lower than ours. 



The only solution we know of to fit the elastic constants 
exactly is to employ a spline form for the pair potential 
as is done in the fitting of bond order potentials ; 24 ' 31 ^ 2 
and we are rather reluctant to make such a departure 
from physical intuition. 



2. fee f-Fe 

Because our models were fitted to the bec Fe lattice 
and elastic constants, it is important to focus on the lower 
part of tableHTlwhich deals with the fee phase of Fe. This 
is 7-Fe which is the base for the austenitic steels and the 
crystal structure adopted by pure Fe above 1185°K4& It 
is well knowniZr— that 7-Fe exists in a high spin ferro- 
magnetic and a low spin (approximately non magnetic) 
modification and we show predictions for both phases in 
table HI1 which we compare with LSDA-GGA calculations 
and experimental observations. It is a mark of transfer- 
ability that both models give a good account of each of 
the two fee phases. Neither model fully captures the large 
and negative C or the softening of C44 of the LSDA-GGA 
in the high spin phase; although they are in better ac- 
cord with experiment than the LSDA-GGA, the proper 
comparison is with the 0°K calculations. The elastic soft- 
ening in 7-Fe is consistent with the measured tempera- 
ture dependence of C in the Invar alloys^ therefore it 
is encouraging that our models are able to describe this 
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TABLE II. Calculated properties using parameters from ta- 
ble U They are compared in the right hand column to either 
experimental values or values calculated using LSDA-GGA, 
the latter written in italics. A proper comparison of the co- 
hesive energy, -E co h, with experiment should take account of 
the spin polarization energy of the free atom which is absent 
in the tight binding limit of infinite separation; this energy 
is as much as^S 0.32 Ry so that the calculated i? C oh should 
amount to 0.63 Ry. Hence the apparent better agreement 
of the orthogonal d model is misleading. Both ferromagnetic 
(FM) and non magnetic (NM) fee Fe is included; we compare 
the experimental data to the FM calculations: the lattice con- 
stant is extrapolated to 0°K>^ the elastic constants are taken 
from phonon dispersion curves^ measured at 1428° K which 
is above the Curie temperature (1043°K) although local mo- 
ments are expected to persist^ LSDA-GGA NM values in 
parentheses refer to the low spin phase. 
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FIG. 4. (color online) Pair potentials in the orthogonal d 
(dotted line) and non orthogonal sd (solid line) tight binding 
models. Note how these are negative at some of the critical 
bond lengths. While this helps to stabilize bec against hep, 
the real benefit is in obtaining a correct elastic constant C' . 
For reference below we also show here the Fe-H pair poten- 
tial in blue (see table [Vj . Vertical lines are placed at first 
and second neighbor distances in bec Fe and at first neighbor 
distances in hep Fe. 
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TABLE III. Calculated surface energies in J/m 2 . Values 
in parentheses are for truncated bulk (unrelaxed) surfaces. 
LSDA-GGA calculations are taken from Spencer et aZ.— 



from data extrapolated from 3°K to 0°K by Adams et aim 
Reference 
Reference 



model 


orthogonal d 


non orthogonal sd 


GGA 


(110) 
(001) 

(111) 


1.77 (1.77) 
2.12 (2.15) 
3.54 (3.85) 


1.53 (1.56) 
1.74 (1.79) 
2.80 (3.34) 


2.27 (2.27) 
2.29 (2.32) 
2.52 (2.62) 



important physical phenomenon at least in principle. It 
has already been shown that elastic and phonon soften- 
ing with increasing temperature in a-Fe is captured in 
the tight binding approximation! 52 1 53 



3. Surface energies 

The proper test of transferability is to carry the mod- 
els into situations of over or under coordination. Here, 
we do this by addressing the surface energies of pure Fe. 
We have set up the (110), (001) and (111) surfaces of 
bec Fe and relaxed the atom positions by energy mini- 



mization using the Hellmann-Feynman forcesj&&££ The 
resulting energies are shown in table Mil in order of de- 
creasing coordination, the most close packed surface be- 
ing (110). We achieve modest, but satisfactory agreement 
with published LSDA-GGA calculations^ at least for the 
two most close packed surfaces. It is in fact notable that 
the LSDA-GGA predicts all the surfaces to have nearly 
the same energy with (111) being a little higher. This is 
not reflected in the tight binding models, indicating lim- 
its to their transferability. The orthogonal d model gives 
the greater spread in energies, demonstrating to some ex- 
tent the greater transferability afforded by the inclusion 
of an s orbital. It is gratifying that both models give 
a qualitative account of surface energies without having 
been fitted, at least in the case of the (110) and (001) the 
latter being of most importance as it's the usual cleavage 
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TABLE IV. Vacancy formation energy, El, in eV, of pure Fe, 
calculated with the orthogonal d and non orthogonal sd tight 
binding models and compared to published LSDA-GGA and 
experimental results. 



model 
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1.61-1.75^ 1.59^ 








2.09^ 


2.0 ± 0.2^ 


unrelaxed 


2.42 


1.36 


2.24^ 2.6(£ 





a Reference 60] 

b Reference 6U 

c Reference 6^ 

d Muon spin rotational 6 ! 

e Quenching-in and electrical resistivity 6 ! 6 ! 

f Positron annihilation, 6 ! but Seeges 6 ! asserts that El S 1.85 eV 

facei 46 ' 56 ' 57 It is also significant in the present context 
that the effect of H on pure Fe and Fe-Si is to enable 
cleavage also on the (110) planesi 5 ! 

4- Vacancy formation energy 



A further test of the transferability is to predict the for- 
mation energy of a vacancy. We do this by constructing 
54 and 53 atom "supercells" of bec Fe (3x3x3 cubic two- 
atom unit cells), one of which has an atom missing. The 
structure is relaxed by energy minimization; its result- 
ing total energy is denoted E(Fesa). The energy of the 
54 atom supercell is denoted E(Fe54). Then the vacancy 
formation energy, neglecting volume relaxation, is 5 ! 

53 

El = E(Fe 53 ) - — £(Fe 54 ). 

Our results are shown in table HVl which also gives values 
for the "unrelaxed" vacancy. As for the surface energies, 
E£ is underestimated by the non orthogonal sd model 
and overestimated by the orthogonal d model. The likely 
error compared to experiment in the latter however is 
more than twice that of the non orthogonal sd model, 
again demonstrating some benefit in transferability of in- 
cluding the non orthogonal s-orbitals. 

IV. ADDING Fe— H INTERACTIONS 

As emphasized before, we will keep the parameters of 
pure Fe unchanged as we seek a model for H in Fe. We 
will find such a model by comparison with properties of 
iron monohydrides of stoichiomctry FcH, that is, the con- 
centrated limit and then test our model's transferability 
into the dilute limit. 

In a series of three papers , 50 ' 68 ' 69 Elsasser et al. have 
made a comprehensive study of the compound FeH in 
the framework of density functional theory. One is in- 
terested in four putative phases, namely fee and bec Fe 



FIG. 5. (color online) To illustrate the tetrahedral (T) and 
octahedral (O) interstices in the bec (upper figure) and fee 
(lower figure) crystals. Note that in the bec lattice the octahe- 
dral site is at the center of a distorted octahedron, unlike the 
fee where it is regular. The distance to the two apical atoms, 
shown here as a horizontal bond, is shorter by a factor l/s/2 
than the distances to the equatorial atoms, two of which are 
shown here in the upper face. This leads in general to the well 
known tetragonal distortion of the bec lattice near octahedral 
interstitial atoms, for example in martensite. For details see 
refs [1^ and [67j • Neither is the tetrahedral interstitial site in 
the bec lattice regular — indeed both octahedral and tetrahe- 
dral bec interstices have tetragonal symmetry. The fee crystal 
structure with all the octahedral sites occupied becomes that 
of cubic rocksalt adopted by many transition metal carbides 
and nitrides. In fee, the tetrahedral site is regular; when half 
these sites are occupied the resulting crystal structure is that 
of zincblende. 




each having one H atom in either tetrahedral or octa- 
hedral sites. These are illustrated in Fig. [5] and Fig. [6] 
shows energy-volume curves for these four phases cal- 
culated using LSDA-GGA in the full potential LMTO 
method^ (see also Fig. 5 in ref (50T|). 

Examination of the upper sketch in Fig. [5] shows that 
the displacement of the tetrahedral interstitial atom in 
the bec structure towards the octahedral site brings the 
impurity atom from above the second neighbor bond, at 
right angles until it finally rests at the bond center. This 
is precisely the situation envisaged by Haas et a/.— in 
their proposal of the screening function, and we therefore 
expect for a model to be transferable, we will require it to 
be non orthogonal. There is also a strong argument for 
the retention of the Fe 4s orbital even though, as we have 
seen, it does not lead to a significantly better model for 
pure Fe than the orthogonal d! 5 - The argument for its in- 
clusion follows from an examination of Fig. [7] which shows 
LSDA-GGA energy bands for bec tetrahedral FeH. The 
bands are colored according to the eigenvector weights 
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FIG. 6. (color online) Cohesive energy and magnetic moment as a function of volume per Fe atom in the four FeH phases 
calculated within the LSDA-GGA (left). Dotted lines denote non magnetic phases. The cohesive energy is with respect to solid 
a-Fe and molecular H2 also calculated using the same energy functional and hence is an approximation to the heat of formation. 
Note that on this basis none of the phases is expected to exist. On the right we show the same quantities calculated in the non 
orthogonal sd tight binding model. We expect that the almost exact degeneracy of bcc TET and fee TET is accidental. 
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coming from LMTO's from H Is (red) or Fe 3d (blue). 
The His band is split off from the Fe 3d bands and has 
similar width. The Fe 4s band which in pure Fe has its 
bottom below the Fe 3d bands and which hybridizes with 
them (see Fig. [T]) is pushed up above the top of the Fe 3d 
bands by repulsion of the H Is band. This means that 
the Fermi energy remains near where it is in pure Fe. 
Roughly speaking one might say that the single 4s elec- 
tron per atom in pure Fe is transferred to the hydrogen 
atom to complete its Is shell, or rather to fill the H Is 
band. At first glance it may seem natural to neglect the 
Fe 4s bands in FeH. But a difficulty will arise if we adapt 
a d-only model by adding just an extra H Is orbital. Hy- 
drogen brings one electron with it and to fill the split-off 
His band an electron will be drawn down from the Fe 3d 
bands consequently lowering the Fermi level. If we were 
only interested in FeH then we could just adjust Nd, the 
number of d electrons; but this will introduce an inconsis- 
tency in going to the dilute limit: will somehow need 
to be continuously adjusted at Fe atoms successively fur- 



ther away from an impurity H. It is very hard to see how 
this problem could be overcome except possibly by al- 
loting two electrons to the hydrogen impurity; while it 
is solved naturally by the Fe 4s falling back into place 
as an Fe atom finds itself remote from the influence of 
impurity. We emphasize that in the non orthogonal sd 
model and its extension to impurities the number of elec- 
trons is not a parameter — as long as all occupied bands 
are included in the hamiltonian we can happily take the 
number of electrons from the periodic table. 

Therefore we take over the pure Fe non orthogonal sd 
model and we add parameters to account for the addi- 
tional H s band. We need Fe-H sda and sscr hopping 
and overlap parameters but we do not require H-H in- 
teraction parameters since even the closest approaching 
interstitial sites are distant more than three times the 
length of the H2 molecular bond. The sda and sscr in- 
tegrals establish the width of the H s band while its po- 
sition with respect to the d bands is set by the on-site 
energy, e s of the H s orbital. We also require Hubbard- 
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FIG. 7. (color online) Energy bands for bcc tetrahedral FeH, calculated at the lattice constant of pure bcc Fe. The upper panels 
show majority and the lower minority spin states. The coloring is such that H-s character is red and Fe-d character is blue. 
Fe-s bands are green. The Fermi energy is indicated by a horizontal line. Note that the Fe-4s band has been pushed above the 
d-bands. Bands on the left are from our tight binding model and on the right are bands calculated in the LSDA-GGA. 
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U parameter o i for H and Fe, but these are not criti- 
cal and 1.2 Ry and 1 Ry are good choices. Essentially 
these lead to approximate charge neutrality as expected 
in metals and their alloys.— For simplicity we take the 
Stoner parameter for H to be zero. Tetrahedral bcc FeH 
is ferrimagnetic, both in LSDA and in our tight binding 
model, the H atom carrying a small moment, less than 
1 /ib ( alig ned opposite to that of the Fe atom cf., Fig. 8 
in ref [If). 

To find the additional parameters we have resorted 
to fitting these to the four equilibrium atomic volumes 
and three cohesive energy differences marked with dashed 
lines on Fig. [5] We do this using Schwefel's multimem- 
bered evolution strategy! 71 1 72 For the Fe-H pair potential 
we employ 




The resulting parameters are displayed in table [V] and 
the hopping integrals are shown graphically in Fig. [8] to 
illustrate their relative magnitudes and ranges. In the 
same figure we show the hopping integrals for Fe which 
are, of course, identical to those of our non orthogonal sd 



TABLE V. H on-site, and Fe-H interaction parameters of 
our tight binding model. All quantities are given in atomic 
Rydberg units. For all these integrals we use ri = 0.8 and 
r c — 2 in units of the pure Fe bcc lattice constant, a = 2.87 A; 
for the pair potential we use n = 0.8 and r c — 0.95 in the 
same units. 
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FIG. 8. (color online) Hopping and overlap integrals as func- 
tions of bond length, r, in the sd non orthogonal model. Ex- 
cept in the case of ddS the dotted lines are the overlap in- 
tegrals corresponding to the hopping integrals of the same 
color. Vertical dotted lines indicate the Fe-H bond length 
in bcc tetrahedral FeH at equilibrium volume, and the Fe-Fe 
bond lengths of the first six neighbors in pure bcc Fe. 
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model of section IIII1 With reference to our remarks in 
section lLTAl we note that all our hopping and overlap inte- 
grals have the simple exponential form up to the distance 
7"i beyond which they are augmented so as to go continu- 
ously and differentiably to zero at r c . These distances arc 
not strictly parameters of the model and are not used in 
the fitting. They are chosen intuitively; for example one 
expects just first neighbors in hep and fee, and first plus 
second neighbors in the bcc structures to be interacting 
through dd hopping whereas the s electrons in pure Fe 
are essentially free electron like and hence "do not take 
kindly to being treated within a TB framework" They 
are best represented by longer ranged interactions. These 
points are illustrated in Fig. [5] and the values of r\ and 
r c can be found in tables HI and fVl The use of fifth degree 
polynomials to augment the tails is necessary to acheive 
a smooth join; it can lead to small kinks as seen in Fig. HI 
but these are designed to fall in between neighbor shells 
and so minimize their effect. This is why the parameters 
ro and r c are made to scale with the lattice constant. The 
resulting energy bands are plotted in Fig. [7] for compar- 
ison with the LSDA-GGA bands. The resulting energy 
volume curves are shown in Fig. [6] The TB model does 
not reproduce the magnetic moments of the LSDA-GGA 
in Fig. [pj quantitatively since this is a sensitive function 
of the density of states at the Fermi level in the non mag- 
netic crystal and our energy bands are only in qualitative 
agreement with the LSDA-GGA. 

Table [VTl summarizes the equilibrium properties of the 
four hydride phases shown in Fig. O The question of site 
selectivity, especially in bcc Fe is important and we will 



TABLE VI. Equilibrium volumes per Fe atom and cohesive 
energies of the four FeH phases following evolution optimi- 
sation, compared to the target values. Cohesive energies are 
relative to the fee octahedral (rocksalt) phase. The final col- 
umn shows the radius of the interstitial site based on a lattice 
of hard spheres at the equilibrium volume of pure Fe and 
taken from Leslie.— All quantities are given in atomic Ryd- 
berg units. 
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revisit it in the dilute limit, below, in section LVB II 



V. PREDICTIONS OF THE Fe-H MODEL 

A. Iron hydride 

Our first test of the tight binding model is to compare 
the resulting adiabatic potential surface section with the 
results of calculations by Elsasser et a/. 69 ' 73 which were 
made in the local density approximation (LDA) to DFT. 
In these calculations the H sublattice is displaced with 
respect to the Fe sublattice in both bcc and fee FeH in 
a chosen set of directions so as to explore the curvatures 
and barriers of the potential energy landscape. For the 
case of the bcc structure, Fig. |pj shows some of the dis- 
placment paths. The potential sections from previous 
LDA 69 and our present tight binding model are shown in 
Fig. [TU] Whereas the relative energies of the tetrahedral 
and octahedral sites have been established by the fitting, 
the remainder of the these curves amount to predictions 
of the tight binding model. They turn out to be be in 
remarkable, quantitative agreement with the LDA calcu- 
lations in the bcc and fee case, the latter being shown 
in Fig. [TTJ These curves exploit to the full the notion 
discussed in section III B[ above, of environment depen- 
dent screening of hopping integrals as the hydrogen ap- 
proaches Fe-Fe first and second neighbor bonds and in- 
deed penetrates the bond to lie directly in between the 
two atoms. It is exactly in this situation that one ex- 
pects the Fe-Fe bond integrals to be strongly modified 
by screening, and clearly our model captures this well in 
a non orthogonal two center description. In particular 
note, in reference to Fig. QJJ that the minimum energy 
(saddle) point along the [101] path lies to the left of the 
point "S" in both LDA and in our TB model. This im- 
plies that the [101]* minimum energy diffusion path in 
reality is bowed slightly towards the center of Fig. [UJ The 
strongest test of the environment dependence however is 
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FIG. 9. (color online) Illustrates the translations of the bcc 
intersitials in constructing our adiabatic potential surfaces, 
after the three dimensional drawing of Fig. 1 in Krimmel et 
aZ.— The figure represents an (010) face of the bcc lattice 
with Fe atoms as black circles at each corner. The octahedral 
sites are shown as squares, the central, filled one being the 
one occupied in octahedral FeH. Of the four tetrahedral sites, 
shown as triangles, one is occupied in tetrahedral FeH and 
this is shown filled in here. The point, S, is midway between 
two tetrahedral sites — the expected diffusion path of H in 
Fe 7 -^ which is highlighted in red here and in Fig. 1101 Those 
displacements which are in the (010) plane are indicated. 



Z 




in the fee hydrides of Fig. [TT] The energy barrier at the 
maximum of the (110) o path, coinciding with the maxi- 
mum of the (001)t path is perfectly rendered by the TB 
model without having been fitted and this corresponds to 
the extreme instance of screening in which the H atom 
becomes positioned at the center of the first neighbor 
Fe-Fe bond (see Fig. 1, ref [z|). 



B. H in Fe — the dilute limit 

We concentrate on three predicted properties of iron 
in this section. First is the dissolution energy^ or zero 
temperature heat of solution of hydrogen in Fe. Included 
in this study is the matter of the site selectivity. Second 
is the binding energy^ or 0°K segregation energy of H to 
the (001) surface of Fe. Third, and of great importance 
to the question of hydrogen embrittlement, is the binding 
of H atoms to a vacancy in Fe. 



TABLE VII. Dissolution energy, in eV, of H in Fe in both 
tetrahedral (TET) and octahedral (OCT) interstices. Present 
results are marked TB, experimental and LSDA-GGA values 
are taken from Jiang and Carter. 76 





TET 


OCT 


TB 
expt. 
GGA 


0.273 
0.296 
0.19 


0.354 
0.32 



1. Dissolution energy 

Following Ramasubramaniam et alX we construct a 
54 atom supercell as we did in section UlI B 41 and whose 
total energy we denoted EiYe^). We then place a hydro- 
gen atom at either a tetrahedral or an octahedral site and 
minimize the total energy by relaxation. The resulting 
total energies are denoted EiYe^H). We do not allow 
the volume to relax. Then the dissolution energy is2i 

E dis = £(Fe 54 H) - £(Fe 54 ) - ^h 2 (4) 

Our model does not contain H-H interactions, but faux 
de mieux we may take En 2 — —4.75 eV from experi- 
ment or from quantum chemistry 20 ' 76 For each of the 
three calculations we employ a mesh ofl2xl2xl2k- 
points and use first order generalized Gaussian integra- 
tion of the Brillouin zone with a width of 2.5 mRy. 77 
Results are shown in table IVIII These are in remarkably 
good quantitative agreement with both observations and 
LSDA-GGA calculations. In particular we predict the 
tetrahedral site to be preferred over the octahedral, as 
is well established^ We may point out here that this is 
not a trivial result: carbon in contrast, while preferring 
the tetrahedral site in the ficticious bec-based carbide, 
transfers to the octahedral site in the dilute limit.— In 
the effective medium theory, upon which the embedded 
atom potentials (EAM) are based, H prefers the octahe- 
dral sitep 7 ^ 



2. H segregation to the ( 001 ) surface of Fe 

Three binding sites of H to the (001) surface of Fe have 
been identified X These are illustrated in Fig.[l2j We have 
constructed supercells of 2 x 2 x 5 cubic two-atom unit 
cells with three layers of vacuum inserted along the long 
axis. The slab contains 40 Fe atoms and the total energy 
of the fully relaxed supercell is denoted -E surf (Fe 4 o). We 
place one H atom at one of the three adsorption sites in 
Fig. [12] and relax the structure by energy minimization. 
Allowing all atoms to relax we denote the total energy 
E s (Fe4o-ff). The associated "adsorption energy" is^ 

E ads = £(Fe 54 H) - E(Fe 54 ) - \E^ 
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FIG. 10. (color online) Adiabatic potential surface sections of bcc FeH: left LDA, 69 right TB. These curves show the energy as 
a function of the displacement of the H sublattice relative to the Fe sublattice. The curves which start at the point "O" refer to 
displacements from the octahedral site phase; a H atom initially at position [§0|] translates in the directions indicated. Along 
the [001] direction it eventually falls into a vacant tetrahedral site (see Fig. |9j. This curve hence represents the transition to 
the tetrahedral-site phase. Translation along [101] takes the H atom to a position midway between two, vacant, tetrahedral 
sites — this point is marked "S". For a H atom initially occupying a tetrahedral site, translation along [101] moves it to an 
adjacent, unoccupied, tetrahedral site, the half-way point being the same point "S". The translation labels are vectors referred 
to Fig. [9] For each case, LDA and TB, the calculations are at fixed atomic volume, namely the equilibrium volume of the bcc 
tetrahedral phase of FeH, see table fVTl ao is the corresponding equilibrium lattice constant. 




and by combining the previous two equations the "bind- 
ing energy" is^ 

^bind = Ea is — -E ac js (5) 

in which the reference energy, or chemical potential, of 
gaseous H 2 has canceled. E^ is is the dissolution en- 
ergy © at a tetrahedral site (table IVIIf) . We have 
calculated the three quantities using a 12 x 12 x 1 k- 
point mesh and the same Brillouin zone integration as 
above. In table I Villi we show our calculated binding en- 
ergies, the displacement 8 in Fig. [12] and the height, h, 
from the (001) surface constructed as the difference in 
z-coordinates of the H atom and the average from the 
four topmost Fe atoms. 



The predictions of our model are only in reasonable 
agreement with the LSDA-GGAi The heights above the 
surface are well rendered; the displacement, 8, is signif- 
icantly larger, but is consistent with the preference for 
tetrahedral site occupancy. As we point out in the cap- 
tion to Fig. [12] 8 = 0.25a = 0.71 A puts the H atom 
into a surface tetrahedral site and our model does exactly 
that; in contrast the LSDA-GGA quite surprisingly re- 
sults in a much smaller 8. In the same vein, the height 



TABLE VIII. Predicted structure and energetics of H ad- 
sorbed on Fe (001). We show for the QT, hollow (H) and 
bridge (B) sites of Fig. [12] the displacent 8 and height, h, 
above the surface (all in A) and the 0°K segregation or bind- 
ing energy, -Ebind, in eV. In parentheses are the LSDA-GGA 
results of Ramasubramaniam et ah 1 





8 


h 


-Ebind 




TB (GGA) 


TB (GGA) 


TB (GGA) 


QT 
H 
B 


0.635 (0.19) 


0.31 (0.38) 
0.27 (0.38) 
0.85 (1.20) 


0.241 (0.768) 
0.191 (0.775) 
0.222 (0.655) 



of the H atom above the bridge site, 0.85 A, is close to 
0.25ao, and we find another local minimum at 0.34 A be- 
low the bridge site. Thus the strongest binding in the 
TB model is to surface tetrahedral sites and the surface 
octahedral site is indeed not a local energy minimum. 
In this way the binding energies are in poor agreement 
with the LSDA-GGA and may reflect the limitations in 
transferability (section fll B[) in that the model retains its 
bulk-like features at the surface. -Ebind is in fact the 0°K 
segregation energy, usually defined as the energy needed 
to remove the impurity from the surface and place it 
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FIG. 11. Adiabatic potential surface sections of fee FeH: left LDA,— right TB. At the point "O" we have the rocksalt phase, 
from which translation of the H sublattice along a (111) direction transforms the structure to the zincblende phase in which 
tetrahedral sites are occupied. The energy maximum between "O" and "T" is located close to where the H atom squeezes 
between an equilateral triangle of Fe atoms in the (111) plane. At the maximum along (110} o , and along (001}t, the H is 
positioned mid- way between two nearest neighbor Fe atoms (see Fig. 1 of ref [zH]). Note that both these two energy barriers 
are predicted by the TB model with quantitative accuracy. The calculations are at the calculated equilibrium volume of the 
fee octahedral (rocksalt) phase of FeH, see table fVTl ao is the corresponding equilibrium lattice constant. 
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FIG. 12. (color online) Three possible binding sites of H on 
the (001) surface of Fe, after Ramasubramaniam et all- Four 
large circles represent the Fe atoms at the corners of a unit 
cell of the (001) face of the bee lattice. At the center is the 
"hollow" site, a smaller circle; this may be displaced along 
[100] by an amount S to become the "quasi-threefold" (QT) 
site indicated by a triangle. The "bridge" site is shown as 
a square. It is important to recognize that the bridge and 
hollow sites in the plane of the truncated bulk surface are 
octahedral interstices, whereas the QT site at S — 0.25ao is 
a tetrahedral site. If a H atom at the bridge site is displaced 
up or down by 0.25ao then it comes to occupy a tetrahedral 
site. Here, ao = 2.87 A is the equilibrium pure a-Fe lattice 
constant. 




• A 




into the interior of the crystal. The LSDA-GGA shows 
the smallest adsorption energy (largest -Ebind) to be at 
the hollow site; whereas we find it at the QT site and 
at this coverage this is not consistent with experiment 
which shows a transition at 100°K from hollow to QT 
site selectivity between about 0.3ML and lML^S while 
our calculations and the LSDA-GGAi are at 0.25ML. 

Both the QT and bridge sites are at local minima in the 
potential energy in our model. This is consistent with the 
LSDA-GGA. 1 However the hollow site is a local saddle 
point having an almost flat energy surface with respect to 
small displacements parallel to the surface; if we displace 
the H atom a sufficient amount then the structure relaxes 
into the QT site occupancy. This is inconsistent with the 
LSDA-GGA in which surprisingly, in view of there being 
another local minimum at QT just 0.19A distant, the 
hollow site is at a local minimum. 1 

To some extent our choice of chemical potential for 
H%, En 2 , is arbitrary; however the observed bond energy 
leads to a very good rendering of the 0°K heat of solu- 
tion (dissolution energy) of H in Fe, table IVIII On the 
other hand it leads to a positive, but small, adsorption 
energy, E^ds, which means that in our model H2 will not 
dissociate on the (001) surface of Fe. In order to model 
the surface adsorption properly we could make an ad hoc 
adjustment of En 2 . This would be at the expense of less 
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accurate Edi S - For example, if we used the Skinner and 
Pcttifor tight binding model of hydrogen^ then we'd 
have Eb_ 2 — —4.30 eV rather than —4.75 eV. In that case 
our dissolution energy in the tetrahedral site becomes 
0.05 eV (rather than 0.27 eV, cf table EH} but the ad- 
sorption energies are then negative as they should be. 
Of course the segregation energies (table IVIIII) remain 
unchanged by this redefinition of the hydrogen chemical 
potential. 



3. H segregation to a vacancy in Fe 

It is believed that the trapping of H to vacancies in 
Fe is of central importance in the effects of H on me- 
chanical behaviori 74 i 81 ' 82 It is also known that dissolved 
hydrogen results in a dramatic increase in the vacancy 
concentration in several metals including Fe ) 83 i 84 caused 
through segregation induced lowering of the vacancy for- 
mation enthalpy 85 We can show that our model is able to 
demonstrate these facts by comparison with LSDA-GGA 
calculations of the 0°K segregation energy, E^ ind (n), of 
up to seven H atoms to a single vacancy in Fe^^ 2 - The 
principal result, which we also predict in our TB model 
is that up to five H atoms may bind to a vacancy with 
a positive segregation energy, but the sixth has a small 
negative E^-An). Here we follow Tateyama and Olmc 62 
and Ramasubramaniam et aim and define E£ id (n) as the 
0°K segregation energy of a H atom from a bulk tetra- 
hedral site to a vacancy to which (n — 1) H atoms are 
already segregated. Hence we set up a 53 atom supercell 
as in section lTlI B 41 then in reference to figure 5 in Rama- 
subramaniam et air- if the vacant site is at [555] in the 

r 1 1 1 



bcc supercell we add H atoms successively in (1) [^t^L 



, 2 ], and (6) [0±±] 



(2) [III], (3) [III], (4) [I i], (5) [U 
octahedral interstices — these are the centers of the six 
{001} faces bounding the vacant site. Finally a seventh 
H atom may be placed at the vacant site. These super- 
cells are relaxed by energy minimization and we denote 
the total energy of the supercell by -E(Fe53H„). Then we 
have^ i 62 ' 86 in analogy with 
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(£(Fe 5 3H, 



-E(Fe53H„„i) — hEn 2 



which is independent of the chemical potential of H. 
Table IIXI shows our segregation energies, compared to 
LSDA-GGA. The relaxation pattern is very simple in all 
cases except n — 3 and n = 5. In the simple instances, 
each H atom relaxes perpendicularly to its {001} face, by 
an amount we denote <5^ ven (n), towards the vacant site. 
The displacement decreases as n increases both in LSDA- 
GGA^ 2 - and our TB model. In each of the cases n = 3 
and n = 5 there is one H atom which follows this trend 
whereas the remaining (n— 1) H atoms are displaced both 
towards the vacancy by i5° dd (n) and, by an amount <5||(n) 
in a direction parallel to the {001} face containing the 
site where the H atom was originally placed, in a (100) 
direction. 



TABLE IX. Segregation of H atoms to a vacancy in Fe. We 
show our model's predicted E^ ind (n) compared to LSDA- 
GGA results,— quoted by Ramasubramaniam et al.,^ in eV. 
Also shown are the displacements of the H atom towards the 
vacancy, and away from the octahedral site in the {001} plane 
in which it was originally placed. In cases of higher symmetry 
the displacement of all H atoms is an amount (n) normal 
to the {001} face and towards the vacant site. In the cases 
n = 3 and n = 5 one atom follows this displacement, while all 
those remaining move both perpendicular to the face — by an 
amount <5x dd (n) — and parallel to the face in a (100) direction 
by an amount 811 (n), rather like the knight's move in chess. A 
displacement 5m = 0.25ao = 0.71 A will take the H atom into 
a tetrahedral site. Displacements are given in A. 
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0.559 
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0.19 
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Table IIXI shows very much better agreement with the 
LSDA-GGA than in the case of surface segregration. 
This probably reflects the better transferability into the 
less undercoordinated environment. Our absolute values 
of E^ ind (n) are no more than 50% underestimated while 
the trends are in perfect accord: we observe the increase 
in segegration energy going from n = 1 to n = 2 imply- 
ing that a H atom segregates more readily to a vacancy 
that has already trapped a H atom. We also see that 
up to five H atoms will segregate exothermally to a va- 
cancy, while the sixth segregates endothermically. The 
displacement patterns in the symmetric cases are con- 
sistent in magnitude with the LSDA-GGA^ 2 and follow 
the trend of decreasing <5<[ von with increasing n. For the 
case n = 1 we obtain 5Y en = 0.25 A which agrees well 
with the LSDA-GGA calculated value of 0.22 A^ 2 . An 
experimental estimate of 0.4 ±0.1 A was obtained for 
deuterium in Fe by ion channeling. 87 Effective medium 
theory for n = 1 results in J^™" = 0.5±0.lAin Fe^ and 
0.46 ± 0.07 A in Nb. 89 The octahedral sites in which the 
H atoms are originally placed correspond to the hollow 
sites at the (001) surface, and as in the surface case the 
atoms relax into the vacuum or vacancy and, symmetry 
permitting, laterally towards the tetrahedral positions. 
The interpretation of Tateyama and Ohno22, that there 
is an electrostatic repulsion between H atoms is uncon- 
vincing to us, since we imagine that this will be screened 
by the electrons in the vacant site. We note that in the 
highly endothermic segregation of a seventh H atom to 
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the vacancy there is still an inward relaxation, at least 
in our model, towards the vacant site, now occupied by 
a H atom. However our -Eb ind (7) is more than five times 
smaller than in LSDA-GGA^ 

We should note, as Kirchheim has pointed outj22, that 
the reduction in enthalpy of the impurity by segregat- 
ing to a defect is entirely equivalent to a reduction in 
the defect's enthalpy of formation. Hence ours and the 
LSDA-GGA binding energies of table IIXI are consistent 
with the observed "superabundant vacancy formation" 
in many metals subject to a high hydrogen fugacit y 83 i 84 
(see Fig. 7, ref HI). 



VI. DISCUSSION AND CONCLUSIONS 

We have described simple and robust tight binding 
models for pure Fe, transferable from bcc into hep and 
fee structures and hence able to describe the common 
phases of Fe, a, 7 and e. Furthermore we have included 
a description of the electronic structure of monohydrides 
and this model has been shown to be transferable into 
the dilute limit of interstitial H impurity in Fe. A simple 
orthogonal d-band model is expected to be most appro- 
priate for the pure transition metals and their alloys 3 -^— 
and indeed the addition of s or p electrons does not usu- 
ally result in better energetics 4^ This is confirmed here 
(table |H]) in the case of bulk elastic constants and struc- 
tural energy difference. The only improvement to bulk 
properties arising from the non orthogonal sd model is 
an improved cohesive energy. Vacancy formation and 
surface energies are somewhat improved in the non or- 
thogonal sd model. 

The focus on transferability is made in section IIVI 
where, while not permitting the parameters of the pure 
Fe model to be adjusted, we seek additional parameters 
to describe Fe-H interactions. We give reasons in sec- 
tion IIVI in addition to the transferability arguments for 
choosing to extend the non orthogonal sd rather than the 
orthogonal d model to the description of hydrogen. There 
are only few additional parameters needed (table [Vj) and 
we emphasize that these were fitted to just seven fidu- 
cial points in the LSDA-GGA energy-volume curves for 
four putative iron hydrides (Fig. [B]). Possibly as a conse- 
quence of our adoption of a non orthogonal model both 
for pure Fe and Fe-H, our resulting model predicts calcu- 
lated adiabatic potential surfaces with quantitative accu- 
racy. It is particularly notable that in these tests H atoms 
are brought perpendicularly towards Fe-Fe bonds to the 
point that the H atom comes between the two host atoms. 
This happens in both bcc and fee hydrides; in the latter 
case a H atom also pushes through the triangle of nearest 
neighbour Fe atoms in the (111) plane and the matching 
to the LDA is excellent ( figs. ITUl and ITT]) . 

Our approach has been to find a model purely by refer- 
ence to the concentrated limit of a stoichiometric mono- 
hydride, FeH, and then to test that model into the di- 
lute limit of H in Fe. Therefore all the results in sec- 



tion [VB] are predictions. In contrast, in constructing a 
classical model Ramasubramaniam et alX needed to put 
all the properties that we describe in section IVB1 into the 
training set for the potential. In consequence, the tight 
binding approach cannot hope to reproduce the quanti- 
tative accuracy that is achieved by a well fitted classical 
model. However dissolution energies, site selectivity and 
vacancy segregation are very well rendered in the model. 
Its most obvious shortcoming is in the prediction of ad- 
sorption energies of H on the (001) surface of Fe. The 
absolute cohesive energy is problematic in LSDA ) 38 i 90 
but even more so in tight binding (see the caption to 
table HIl . Possibly for this reason we find that H2 will 
not dissociate on the (001) surface if we use the known 
binding energy of the H2 molecule as our reference. In 
future work we will need to account for molecular hydro- 
gen and this matter will be revisited. On the other hand 
qualitatively the TB model gives a reasonable account 
of H adsorption which is certainly a subtle and complex 
problem in surface physics. In this way the TB model 
does not transfer faultlessly into the problem of surface 
energetics. Our predictions of segregation to a vacancy, 
in contrast, are in very good accord with the known the- 
oretical LSDA-GGA results and experimental facts. In 
particular we predict that a vacancy will bind up to five 
H atoms exothermically and that the segregation energy 
is somewhat larger to a vacancy at which one H atom 
is already bound. The trapping of vacancies is central 
to the mechanism of the action of H on the mechanical 
properties of Fe alloys / 4 ' 81 ' 82 

In conclusion, the quantum mechanical tight binding 
approximation lies between the first principles LSDA and 
the atomistic classical approach to defect energetics in 
iron. Because the TB approximation is grounded in elec- 
tronic structure theory it may be applied to this ques- 
tion rather easily and just a few parameters — adjustable 
within intuitive limits — are required. Because of this and 
because of its simplicity the TB approach may give rise to 
a better understanding than the LSDA, which after much 
labor produces a total energy and force, often without 
clear insight to their origins. In contrast the huge num- 
ber of parameters and the rather opaque functional form 
of the interatomic interactions in the classical potentials, 
while able to model many properties quantitatively, must 
be at risk of failure once they are transferred into situa- 
tions for which they were not fitted. Therefore we expect 
the TB approximation to provide a useful and comple- 
mentary tool to the classical potentials, and once aug- 
mented with parameters to describe carbon, to become 
competitive in the atomistic simulation of the properties 
of iron and steel. 
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